

/*******************************
********************************
Replication file
Lahore Bus Rapid Transit analysis
********************************
*******************************/

********************************
********************************
* Set globals 
********************************
********************************

	global replication ""X:\Dropbox\Metrobus study\Replication""
	global repout ""X:\Dropbox\Metrobus study\Replication\tables\""
	
	global covariates gender ageyears ageyearssq yearsed yearsedsq dist_cent dist_cent_sq census_* 
	

********************************
********************************
* HH data analysis
********************************
********************************

	cd $replication

	use metrobus_panel_repfile, clear

	
	****************
	* Summary stats
	****************
			
	cd $repout
	estpost summarize workoutside ///
		vehiclecommute ///
		publiccommute publiccommute_cond  ///
		gender ageyears yearsed dist_cent ///
		if year==2015

	esttab using sumstats.csv, ///
		cells("count mean(fmt(3)) sd(fmt(2)) min max") ///
		label csv replace ///
		title("Summary statistics - HH survey") 
	
	
	****************
	*Check for baseline balance 
	****************

	global balance "" 
	
		foreach i in publiccommute avgtime vehiclecommute motorcycles {
		
		reg `i' treat1 dist_anystop diststopgroup_* ///
			$covariates [pw = weightperhh] ///
			if year == 2012 , cluster(zone_id) 
			est sto `i'2
			sum `i' if year==2012 & treat1==0 & treat2==0
			estadd scalar cmean = round(r(mean), .01)
			estadd local donut "Yes"
			estadd local covariates "Yes"
			
		reg `i' treat1 treat2 dist_anystop diststopgroup_* $covariates ///
			[pw = weightperhh] if year == 2012 , ///
			cluster(zone_id) 
			est sto `i'3
			sum `i' if year==2012 & treat1==0 & treat2==0
			estadd scalar cmean = round(r(mean), .01)
			estadd local donut "Yes"
			estadd local covariates "Yes"
			
		
		global balance $balance `i'2 `i'3 

		}

		cd $repout 
		esttab $balance using balance2012.csv, ///
			replace csv ///
			star(* .1 ** .05 *** .01) label se t(4) b(4)  ///
			keep(treat1 treat2) ///
			title("Baseline balance tests - 2012") ///
			scalars("cmean Control group mean" "donut Donut FE" "covariates Additional covariates") 
		
		

	****************
	* First stage
	****************
		
		foreach i in traveltime_ travelfare_ {

		reg `i' mindist1 dist_anystop if post==1 & uniquespn==1, cluster(zone_id)

			est sto first1`i'
			estadd local donut "No"
			estadd local est "XS"

		reg `i' mindist1 mindist1_sq dist_anystop dist_anystop_sq ///
			if post==1 & uniquespn==1, cluster(zone_id)

			est sto first3`i'
			estadd local donut "No"
			estadd local est "XS"

		reg `i' mindist1 mindist1_sq dist_anystop dist_anystop_sq post ///
			mindist1_post mindist1_sq_post dist_anystop_post dist_anystop_sq_post ///
			if uniquespn==1, cluster(zone_id)

			est sto first5`i'
			estadd local donut "No"
			estadd local est "Panel"
			
		reg `i' mindist1 mindist1_sq dist_anystop dist_anystop_sq post ///
			mindist1_post mindist1_sq_post  ///
			dist_anystop_post dist_anystop_sq_post diststopgroup* ///
			if uniquespn==1, cluster(zone_id)

			est sto first6`i'
			estadd local donut "Yes"
			estadd local est "Panel"

		}

		cd $repout

		esttab first1traveltime_ first3traveltime_ first5traveltime_  ///
			first6traveltime_ using firststage.csv, replace csv ///
			star(* .1 ** .05 *** .01) se t(4) b(4) ///
			label drop(_cons diststopgroup*) ///
			nomtitles nodepvars  scalars("donut Donut FE" "est Specification")

		esttab first1travelfare_ first3travelfare_ first5travelfare_  ///
			first6traveltime_ using fare.csv, replace csv ///
			star(* .1 ** .05 *** .01) se t(4) b(4) ///
			label drop(_cons diststopgroup*) ///
			nomtitles nodepvars  scalars("donut Donut FE" "est Specification")

			
			
		ivreg2 commutemode1 mindist1 dist_anystop ///
			[pw = weightperhh] if year==2015, cluster(zone_id)
		est sto takeup0
		estadd local covariates "No" 
		estadd local donut "No" 
	
		ivreg2 commutemode1 mindist1 dist_anystop $covariates ///
			[pw = weightperhh] if year==2015, cluster(zone_id)
		est sto takeup1
		estadd local covariates "Yes" 
		estadd local donut "No" 
		
		ivreg2 commutemode1 mindist1 mindist1_sq dist_anystop dist_anystop_sq ///
			$covariates [pw = weightperhh] if year==2015, cluster(zone_id)
		est sto takeup2
		estadd local covariates "Yes" 
		estadd local donut "No" 

		ivreg2 commutemode1 mindist1 mindist1_sq dist_anystop dist_anystop_sq ///
			diststopgroup_* $covariates [pw = weightperhh] if year==2015, ///
			cluster(zone_id) partial(diststopgroup_*) 
		est sto takeup3
		estadd local covariates "Yes" 
		estadd local donut "Yes" 

		sum commutemode1 if year==2015
		estadd scalar mean = round(r(mean), .01)

		cd $repout		
		esttab takeup0 takeup1 takeup2 takeup3 using takeup.csv, ///
			replace csv ///
			star(* .1 ** .05 *** .01) se t(4) b(4) ///
			title("Takeup" "Dependent variable: Uses green line for daily commute") ///
			drop($covariates) label nomtitles ///
			scalars("mean Sample mean dependent variable" ///
				"covariates Additional control variables" "donut Donut ring FE") 			
			
			
		
	****************
	*Main results 
	****************

			global spec0 dist_anystop dist_anystop_sq ///
				(traveltime_ = mindist1 mindist1_sq) [pweight = weightperhh] ///
				if post==1, ///
				cluster(zone_id) savefirst 
			
			global sample0 post==1 
			
			global spec1 $covariates diststopgroup_* dist_anystop dist_anystop_sq ///
				(traveltime_ = mindist1 mindist1_sq) [pweight = weightperhh] ///
				if post==1, ///
				cluster(zone_id) savefirst partial(diststopgroup_* $covariates) 
				
			global sample1 post==1
				
			
			global spec2 $covariates diststopgroup_* dist_anystop ///
				dist_anystop_sq dist_anystop_post dist_anystop_sq_post ///
				mindist1 mindist1_sq post ///
				(traveltime_ = mindist1_post mindist1_sq_post) [pweight = weightperhh] ///
				if year > 2009, ///
				cluster(zone_id) savefirst partial(diststopgroup_* $covariates) 

			global sample2 year > 2009 
				
				
			global spec3 $covariates diststopgroup_* dist_anystop dist_anystop_sq ///
				dist_anystop_post dist_anystop_sq_post ///
				mindist1 mindist1_sq post ///
				(traveltime_ = mindist1_post mindist1_sq_post) ///
				[pweight = weightperhh] ///
				if year > 2009  & (treat1==1 | treat2==1) , ///
				cluster(zone_id)  savefirst partial(diststopgroup_* $covariates) 
			
			global sample3 year > 2009 & (treat1==1 | treat2==1) 

			
			global spec4 $covariates diststopgroup_* dist_anystop dist_anystop_sq ///
				dist_anystop_post dist_anystop_sq_post ///
				mindist1 mindist1_sq post ///
				(traveltime_ = mindist1_post mindist1_sq_post) ///
				[pweight = weightperhh] ///
				if year > 2009 & (treat1==1 | zonegroup_ref==0) ///
				, ///
				cluster(zone_id)  savefirst partial(diststopgroup_* $covariates) 
			
			global sample4 year > 2009 & (treat1==1 | zonegroup_ref==0) 
			
			global spec5 $covariates diststopgroup_* ///
				dist_anystop dist_anystop_sq dist_anystop_post dist_anystop_sq_post ///
				mindist2 mindist2_sq post ///
				(predict_traveltime_ = mindist2_post mindist2_sq_post) ///
				[pweight = weightperhh] ///
				if year > 2009 & (zonegroup_ref==0 | zonegroup_ref==2), ///
				cluster(zone_id)  savefirst partial(diststopgroup_* $covariates) 

			global sample5 year > 2009 & (zonegroup_ref==0 | zonegroup_ref==2) 
	
	
 
			foreach i in publiccommute vehiclecommute workoutside motorcycles {

			forvalues k = 0 / 5 { 
			
			ivreg2 `i' ${spec`k'} 
				
				est sto `i'`k' 
				estadd local covariates "No" 
				estadd local donut "No" 
				sum `i' if ${sample`k'} 
				estadd scalar premean = round(r(mean), .01)

			} 
			
			cd $repout
			esttab `i'0 `i'1 `i'2 `i'3 `i'4  ///
				using mainestimates_`i'.csv, replace fragment csv ///
				star(* .1 ** .05 *** .01) se t(4) b(4) label keep(traveltime_) ///
				nomtitles nonumbers nolines ///
				scalars("covariates Additional control variables" "donut Donut FE" ///
				"sample Geographic sample" "est Specification" ///
				"jp Hansen's J p-value" "rkf First-stage F-stat" ///
				"premean Sample mean pre") 

			} 
		
	
	

	
	****************
	* Heterogeneity
	****************
			
		global outcomeshetsubset publiccommute
		
		global groups educated 
		
		foreach j in $groups {
		
			foreach i in $outcomeshetsubset {

				forvalues k = 0 / 1 {

				ivreg2 `i' dist_anystop dist_anystop_sq ///
					(traveltime_ = mindist1 mindist1_sq) [pweight = weightperhh] ///
					if post==1 & `j'==`k', ///
					cluster(zone_id) savefirst
				
					est sto `i'0`k'
					estadd local covariates "No" 
					estadd local donut "No" 
					estadd local sample "Full"
					estadd local est "XS"
					sum `i' if year==2012 
					estadd scalar premean = round(r(mean), .01)
				
				ivreg2 `i' $covariates diststopgroup_* dist_anystop dist_anystop_sq ///
					(traveltime_ = mindist1 mindist1_sq) [pweight = weightperhh] ///
					if post==1 & `j'==`k', ///
					cluster(zone_id) savefirst partial(diststopgroup_* $covariates) 
				
					est sto `i'1`k'
					estadd local covariates "Yes" 
					estadd local donut "Yes" 
					estadd local sample "Full"
					estadd local est "XS"
					sum `i' if year==2012 
					estadd scalar premean = round(r(mean), .01)

				ivreg2 `i' $covariates diststopgroup_* dist_anystop ///
					dist_anystop_sq dist_anystop_post dist_anystop_sq_post ///
					mindist1 mindist1_sq ///
					(traveltime_ = mindist1_post mindist1_sq_post) [pweight = weightperhh] ///
					if year > 2009 & `j'==`k', ///
					cluster(zone_id) savefirst partial(diststopgroup_* $covariates) 

					est sto `i'2`k'
					estadd local covariates "Yes" 
					estadd local donut "Yes" 
					estadd local sample "Full"
					estadd local est "Panel"
					sum `i' if year==2012 
					estadd scalar premean = round(r(mean), .01)

				ivreg2 `i' $covariates diststopgroup_* dist_anystop dist_anystop_sq dist_anystop_post dist_anystop_sq_post ///
					mindist1 mindist1_sq (traveltime_ = mindist1_post mindist1_sq_post) [pweight = weightperhh] ///
					if year > 2009 & (treat1==1 | treat2==1) & `j'==`k', ///
					cluster(zone_id)  savefirst partial(diststopgroup_* $covariates) 

					est sto `i'3`k'
					estadd local covariates "Yes" 
					estadd local donut "Yes" 
					estadd local sample "T1 T2"
					estadd local est "Panel"
					sum `i' if year==2012 & (zonegroup_ref==1 | zonegroup_ref==2)
					estadd scalar premean = round(r(mean), .01)

				ivreg2 `i' $covariates diststopgroup_* dist_anystop dist_anystop_sq dist_anystop_post dist_anystop_sq_post ///
					mindist1 mindist1_sq (traveltime_ = mindist1_post mindist1_sq_post) [pweight = weightperhh] ///
					if year > 2009 & (treat1==1 | zonegroup_ref==0) & `j'==`k', ///
					cluster(zone_id)  savefirst partial(diststopgroup_* $covariates) 

					est sto `i'4`k'
					estadd local covariates "Yes" 
					estadd local donut "Yes" 
					estadd local sample "T1 C"
					estadd local est "Panel"
					sum `i' if year==2012 & (zonegroup_ref==1 | zonegroup_ref==0)
					estadd scalar premean = round(r(mean), .01)

				}
			
				cd $repout 
				esttab `i'00 `i'10 `i'20 `i'30 `i'40 using het_`j'_`i'_0.csv, fragment replace csv ///
					star(* .1 ** .05 *** .01) se t(4) b(4) label keep(traveltime_) ///
					title("IV estimates - by `j'") nomtitles ///
					scalars("group Subsample" "donut Donut FE" "sample Gegraphic sample" "est Specification" "jp Hansen's J p-value" "rkf First-stage F-stat"  "controlmean Control group mean") ///
					note("Excludes HHs and individuals who moved in 3 years ago or less.  Standard errors clustered by zone (48 zones).") 

				esttab `i'01 `i'11 `i'21 `i'31 `i'41 using het_`j'_`i'_1.csv, fragment replace csv ///
					star(* .1 ** .05 *** .01) se t(4) b(4) label keep(traveltime_) ///
					title("IV estimates - by `j'") nomtitles ///
					scalars("group Subsample" "donut Donut FE" "sample Geographic sample" "est Specification" "jp Hansen's J p-value" "rkf First-stage F-stat"  "controlmean Control group mean") ///
					note("Excludes HHs and individuals who moved in 3 years ago or less.  Standard errors clustered by zone (48 zones).") 

				}
				
				
		}

	estimates drop *
	
		
	****************
	*Binary treatment variable 
	****************
		ivreg2 	publiccommute treat1 treat1_post diststopgroup_* ///
				post ///
				[pweight = weightperhh] ///
				if year > 2009 , ///
				cluster(zone_id) 
		
		est		sto binary1
		estadd 	local subsample "All" 
		estadd 	local covariates "No" 
		estadd 	local donut "Yes"
		estadd	local sample "T1 T2 C" 
		sum 	publiccommute if year==2012 			
		estadd 	scalar premean = round(r(mean), .01)
						
		ivreg2 	publiccommute treat1 treat1_post $covariates diststopgroup_* dist_anystop ///
				dist_anystop_sq dist_anystop_post dist_anystop_sq_post ///
				post ///
				[pweight = weightperhh] ///
				if year > 2009 , ///
				cluster(zone_id) 

		est		sto binary2
		estadd 	local subsample "All" 
		estadd 	local covariates "Yes" 
		estadd 	local donut "Yes"
		estadd	local sample "T1 T2 C" 
		sum 	publiccommute if year==2012 			
		estadd 	scalar premean = round(r(mean), .01)
		
		cd 		$repout
		
		esttab	binary1 binary2 using binaryspecification.csv, replace ///
				star(* .1 ** .05 *** .01) se t(4) b(4) label keep(treat1 treat1_post) ///
				nomtitles nodepvars csv ///
				scalars("covariates Additional control variables" ///
				"donut Donut FE" "sample Geographic sample" ///
				"premean Sample mean pre") ///
			
	
		
		
****************
* Community panel 
****************
		
	cd $replication 

	use community_panel_repfile, clear
		
	cd $repout
	
	twoway (kdensity traveltime_ if year==2012 & mindist1 < 5) ///
		(kdensity traveltime_ if year==2015 & mindist1 < 5), ///
		legend(label(1 "Pre BRT") label(2 "Post BRT")) ///
		title("Pre and post BRT public transport access" "Time to reach center by public transport" "Treatment areas") ///
		note("Notes: Public transport access variable is community respondent's report" "of fastest travel time to central point (Kalma Chowk)" "using only BRT, bus, wagon and/or walking." "Sample is all areas < 5km from a BRT stop.")
		graph export traveltimedisttreat.png, replace 
	
	twoway (kdensity traveltime_ if year==2012 & mindist1 >= 5) ///
		(kdensity traveltime_ if year==2015 & mindist1 >= 5), ///
		legend(label(1 "Pre BRT") label(2 "Post BRT")) ///
		title("Pre and post BRT public transport access" "Time to reach center by public transport" "Control areas") ///
		note("Notes: Public transport access variable is community respondent's report" "of fastest travel time to central point (Kalma Chowk)" "using only BRT, bus, wagon and/or walking." "Sample is all areas > 5km from a BRT stop, including T2 (planned) and C (cancelled) lines.")
		graph export traveltimedistcontrol.png, replace 

	twoway (kdensity travelfare_ if year==2012 & mindist1 < 5) ///
		(kdensity travelfare_ if year==2015 & mindist1 < 5), ///
		legend(label(1 "Pre BRT") label(2 "Post BRT")) ///
		title("Pre and post BRT public transport access" "Fare to reach center by public transport" "Treatment areas") ///
		note("Notes: Public transport access variable is community respondent's report" "of fare in PKR for fastest trip to central point (Kalma Chowk)" "using only BRT, bus, wagon and/or walking." "Fare adjusted to 2015 PKR." "Sample is all areas < 5km from a BRT stop.")
		graph export travelfaredisttreat.png, replace 
	
	twoway (kdensity travelfare_ if year==2012 & mindist1 >= 5) ///
		(kdensity travelfare_ if year==2015 & mindist1 >= 5), ///
		legend(label(1 "Pre BRT") label(2 "Post BRT")) ///
		title("Pre and post BRT public transport access" "Fare to reach center by public transport" "Control areas") ///
		note("Notes: Public transport access variable is community respondent's report" "of fare in PKR for fastest trip to central point (Kalma Chowk)" "using only BRT, bus, wagon and/or walking." "Fare adjusted to 2015 PKR." "Sample is all areas > 5km from a BRT stop, including T2 (planned) and C (cancelled) lines.")
		graph export travelfaredistcontrol.png, replace 
		
	
****************
* Descriptives 
****************
		
	cd $replication

	use ridersurvey_repfile , clear 


	cd $repout 

	graph bar pre_cost_real post_cost [aweight = weightresponders], ///
		title("BRT rider survey: Total fare for the same Green Line trip" "before and after Green Line") ///
		ytitle("Total fare in 2015 PKR; 100 PKR = 1 USD") ///
		legend(label(1 "Trip cost on other modes" "before Green Line (2012)") ///
		label(2 "Total fare today" "including travel" "to and from Green Line")) ///
		note("Notes: Difference is significant at the 0.01% level.")
		graph export pre_post_cost.png, replace

	graph bar [aweight = weightresponders], over(Q5) title("BRT Rider survey:" "How did you get to the station?")
		graph export firstmode.png, replace 

	graph bar if Q5==1 [aweight = weightresponders], over(timecat) ///
		title("How long did you travel to the station?" "Respondents who traveled to station by foot")
		graph export timetostationfoot.png, replace 
	
	graph bar if Q5 !=1 [aweight = weightresponders], over(timecat) ///
		title("How long did you travel to the station?" "Respondents who traveled to station by vehicle")
		graph export timetostationmotor.png, replace 
	
	graph bar foot rickshaw qingqi buswagon cycle motorbike car if none==0 [aweight = weightresponders], ///
		percent title("BRT Rider survey:" "How did you make this trip previously?") ///
		legend(label(1 "Walking") label(2 "Rickshaw") label(3 "Qingqi") label(4 "Bus / wagon") label(5 "Cycle") label(6 "Motorbike") label(7 "Car"))
		graph export previousmode.png, replace 

	graph bar public private both if none==0  [aweight = weightresponders], ///
		percent title("BRT Rider survey:" "How did you make this trip previously?") ///
		legend(label(1 "Public transport only") label(2 "Private transport only") label(3 "Both"))
		graph export previousmode2.png, replace 
	
	graph bar public private both if none==0  [aweight = weightresponders], over(ed) percent stack ///
		title("BRT Rider survey:" "How did you make this trip previously?") ///
		legend(label(1 "Public transport only") label(2 "Private transport only") label(3 "Both"))
		graph export previousmode_ed.png, replace 


	reshape long wtp_yes_, i(id) j(fare) 

	gen wtp_yes_100 = wtp_yes_ * 100

	graph bar wtp_yes_100 [aweight = weightresponders], over(fare) ///
		title("BRT Rider survey:" "Hypothetical Willingness to Pay Higher Fares" "Fares in PKR (105 PKR = 1 USD)") ///
		ytitle("Percentage of riders who say they are willing to pay" "hypothetical increase in fare for same trip") 
	graph export wtp.png, replace 
		
		
		
	cd $replication

	use his_treatment_trip_repfile , clear 

	cd $repout 
		
	graph bar public private walk , over(educ_level_recode, label) percent stack ///
		title("Commuting on public and private transport" "by education level pre BRT" "2010 HIS survey")  ///
		legend(label(1 "Public") label(2 "Private") label(3 "Walking") )
		graph export modes_ed_his.png, replace 
		
	graph bar mode_trip_simple_* , over(educ_level_recode) ///
		percent stack title("Commuting on public and private transport" "by education level pre BRT" "2010 HIS survey") ///
		legend(label(1 "Walking") label(2 "Bicycle") label(3 "Motorcycle") label(4 "Car") label(5 "Public bus / wagon") label(6 "Rickshaw / motorcycle rickshaw") label(7 "Other")) 
		graph export modes_ed_his_detail.png, replace 
		
		
	cd $replication

	use his_rider_appended_repfile , clear 

	cd $repout 
		
	collapse edu_comparable_* [aweight = weightresponders], by(source) 
	
	graph bar edu_comparable_*, over(source) percent stack ///
		legend(label(1 "Less than HS") label(2 "HS degree (FA / FSC)") label(3 "University") label(4 "Postgraduate")) ///
		title("Education levels of BRT riders" "and pre-BRT large bus riders in same areas" "Work commute only") ///
		bar(1, color(gs10)) bar(2, color(gs6)) bar(3, color(gs4)) bar(4, color(black)) 
	graph export edu_comparison_workcommute.png, replace 
